home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Shareware Grab Bag
/
Shareware Grab Bag.iso
/
090
/
byt1186b.arc
/
RUNGKUT.LBR
/
RUNKUT.DOC
< prev
next >
Wrap
Text File
|
1986-04-11
|
29KB
|
794 lines
1. RUNKUT - AN INITIAL VALUE ORDINARY DIFFERENTIAL EQUATION
SOLVER.
1.1. INTRODUCTION
RUNKUT is a subroutine that contains three Runge-Kutta initial
value ordinary differential equation (IVP ODE) solvers They are:
- A fourth order Fehlberg method.
- A fifth order Verner method.
- An seventh order Verner method.
RUNKUT solves coupled IV ODE's of the form:
Y' = F(x,Y)
There is no limit to the number of coupled equations in a system
to be solved other than the amount of memory available. The
differential equations are defined in a user-supplied
subroutine--more about that later. The subroutine is called
using the following FORTRAN call statement:
CALL RUNKUT(XA,Y,XB,NEQN,TOLA,TOLR,HSTART,WORK,
& IMETH,IERROR,ICOM,FUNC)
and must be called from a main "calling" program.
1.2. VARIABLES PASSED TO RUNKUT.
The variables passed are explained below and are listed in the
order in which they appear in the call statement. The
superscripted numbers refer to notes appearing at the end of this
section.
XA (Real-Input)
XA is the starting point of the interval over which
integration of the ODE's is to be performed.(1)
Y (Real array of size NEQN - Input/Output)
On entry, Y must contain the initial values of the
Y-variables--that is Y must contain the solutions at the
point XA. On exit from RUNKUT, Y will contain the solutions
at the end of the specified interval--that is at XB.(2)
XB (Real-Input)
XB specifies the end of the interval over which integration
is to be performed--the point at which solutions to the
differential equations are desired.(1,2)
NEQN (Integer-Input)
NEQN specifies the number of coupled differential equations
in the system to be solved.
TOLA (Real-Input)
TOLA is the absolute error tolerance required on the
solution.(3)
TOLR (Real-Input)
TOLR is the relative error tolerance required on the
solution.(3)
HSTART (Real-Input/Output)
On the very first call to RUNKUT in a calling sequence,
HSTART must contain a non-zero value as an initial estimate
for the step-size.(1) The actual value used is relatively
unimportant as the program will optimize step-sizes to keep
the tolerances within specified limits and keep the number of
function evaluations to a minimum. It is left to the user's
discretion to choose a value that will not require excessive
readjustment by the program. On exit from RUNKUT, this
variable will contain the last step size used by the program.
It is recommended that for subsequent calls to RUNKUT, the
step size returned by the program be used as the value for
HSTART to solve the equations over the immediately following
interval.
WORK (Real array of minimum size NEQN*17 - Input)
This is a work area made available to RUNKUT by the calling
program.
IMETH (Integer-Input)
Indicates to the solver, the order of the method to be used.
IMETH can take on the following values:
1. - Fourth order Fehlberg method
2. - fifth order Verner method
3. - seventh order Verner method
While solving a given set of differential equations it is
possible to change the order of the method over sucessive
intervals simply by setting IMETH to the desired value and
resetting ICOM(1) to 0 on each subsequent call to RUNKUT in
which a change of order is required.(4)
IERROR (Integer-Output)
IERROR is an error status indicator returned by RUNKUT after
each call. It can have one of the following values:
0 - no error.
1 - the sum of TOLA and TOLR is less than 100 times machine
epsilon. This can be caused by setting both TOLA and TOLR
to 0, or by specifying values that are too small for the
program to handle sucessfully. In this case TOLR is set
to a default value.
2 - the problem is too stiff for the method to handle, or
there is a discontinuity in the function. It is
impossible to proceed in either case.(2)
3 - both conditions 1 and 2 above have occured.
IERROR is set to 0 on every entry to RUNKUT.
ICOM (Integer array of size 4)
ICOM is a communications vector. Each array element is used
to pass a specific item of information (as defined below) to
and from the subroutine.
ICOM(1) (Input)
ICOM(1) must be set to 0 on the first call to RUNKUT.
ICOM(1) is internally set to a non-zero value by the program
after the first call. For subsequent calls to RUNKUT with the
same set of equations, ICOM(1) must be left at this non-zero
value. It is possible to change the order of the method used
(see also IMETH) over sucessive intervals. In this case
ICOM(1) must be reset to 0 each time the order of the method
is changed.(4)
ICOM(2) (Input)
A flag that indicates to the program whether to perform
checks on the specified values of TOLA and TOLR.
0 - no checking of the error tolerances.
1 - program checks if the sum of TOLA and TOLR is at least
100 times machine epsilon. If not, TOLR is set to a
default value and the error indicator IERROR is set to 1.
If ICOM(2) is set to 0, the program assumes that it will be
receiving non-zero values of at least one of TOLA or TOLR.
Failure to check this before entry to RUNKUT could result in
fatal program errors. Also, as will be shown in an example,
decreasing the error tolerances does not always result in
improved solution accuracy.
ICOM(3) (Input)
A flag that indicates to the program whether to use the
default values of mimimum and maximum step size set by the
subroutine.(5)
0 - default values of the maximum and minimum step size are
used.
1 - allows the user to set values for the maximum and minimum
permissible step size. These values are set by including
the following FORTRAN common block statement in the main
"calling" program:
COMMON /CONS/HMIN,HMAX
ICOM(4) (Output)
Status flag that indicates the presence of round-off
error in the solution.
0 - No round-off error.
1 - Severe round-off error possible in answer.
ICOM(4) is reset to 0 on entry to RUNKUT.
FUNC (Function name)
FUNC is replaced by the name of the subroutine that is
part of the main user routine in which the differential
equations are specified (See section on user supplied
routines). The subroutine name must be declared external
by the following FORTRAN statement:
EXTERNAL [SUBROUTINE NAME].
1.2.1. NOTES
1. XA could be greater than XB. In other words it is
possible to use RUNKUT to solve the differential
equations in a negative X-direction given initial values
at XA. In that case the step-size will also be negative.
HSTART could be specified as a negative value, but that
is not essential as the program internally adjusts the
algebraic sign of the step-size.
2. If a discontinuity or extreme stiffness is encountered,
RUNKUT will return in XB, the value of X closest to the
point of discontinuity or stiffness, at which the
solutions are still within the specified error bounds.
The array Y will then contain solutions at this value of
X.
3. If the sum of TOLA and TOLR is less than 100 times
machine epsilon (for example if neither were defined in
the main program), TOLR is set to a default value of
10(6) times machine epsilon. TOLA is not set to any
default value. The program uses a combination of TOLA and
TOLR to determine the accuracies of the solutions, and
specifying TOLR to a certain value almost always gives
better results than specifying TOLA to the same value.
4. Changing the order of the method on every subsequent call
to RUNKUT with a given system of equations can lead to
excessive computaion times since a large number of
constants are re-evaluated each time the order is
changed.
5. The default values used by the subroutine are 10(-8) for
HMIN, and for HMAX as follows:
0.5 - if the fourth order method is used.
1.0 - if the fifth order method is used.
2.0 - if the seventh order method is used.
However, if any of these values are larger than the
absolute value of the interval width, HMAX is set to the
interval width. It is sometimes necessary to specify HMAX
to a small value to keep the step sizes to relatively
small values that keep the runge-kutta methods within
regions of stability. The subroutine detects the onset of
problem "stiffness" or the occurence of a discontinuity
in the function by checking if the step size is smaller
than HMIN. If the equations are being solved in the
negative X-direction, HMIN and HMAX are the smallest and
largest permissible step sizes in the negative
X-direction respectively. In other words, HMIN and HMAX
always define absolute constraints on the step size.
6. Machine epsilon on the IBM PC is very much a function of
the compiler used (even with an 8087 math coprocessor
installed). For example Microsoft's Fortran-77 compiler
(with 8087 support) generates an epsilon of about
10(-16), whereas IBM/Ryan-McFarland's Professional
Fortran-77 generates an epsilon of about 10(-20)
1.3. THE USER-SUPPLIED ROUTINES
Two routines must be supplied:
1. A main "calling" program that calls RUNKUT. It must also
define the initial values of Y at XA, set TOLA and TOLR,
set ICOM(1) to 0 on the first call to RUNKUT, define the
number of equations in the sytstem to be solved,
dimension a WORK array for RUNKUT, and check for the
error status as passed back from RUNKUT through IERROR.
See the examples at the end of this section.
2. A subroutine containing the system of differential
equations. The subroutine name must be declared EXTERNAL
in the main calling program and its name passed to RUNKUT
through the parameter FUNC-see above. The subroutine is
defined using the follwing FORTRAN statement:
SUBROUTINE [FUNC] (X,Y,YPRIME,NEQN)
where the parameters have the following definitions:
X (Real)
Contains the current value of X as set in RUNKUT. Do
not alter this value within [func]
Y (Real - array of size NEQN)
Contains the current solutions at X. Do not alter in
[func].
YPRIME (Real - array of size NEQN)
YPRIME is set in the subroutine to the differentials
dy/dx as follows:
YPRIME(1) = F/1/(X,Y(1),...,Y(NEQN))
.
.
YPRIME(NEQN)=F/NEQN/(X,Y(1),...,Y(NEQN))
where F/i/ is the i'th differential function.
1. RUNKUT - AN INITIAL VALUE ORDINARY DIFFERENTIAL
EQUATION SOLVER.
1.1. INTRODUCTION
RUNKUT is a subroutine that contains three
Runge-Kutta initial value ordinary differential
equation (IVP ODE) solvers They are:
- A fourth order Fehlberg method. -
A fifth order Verner method. - An seventh order
Verner method.
RUNKUT solves coupled IV ODE's of the form:
Y' = F(x,Y)
There is no limit to the number of coupled
equations in a system to be solved other than the
amount of memory available. The differential
equations are defined in a user-supplied
subroutine--more about that later. The subroutine is called
using the following FORTRAN call statement:
CALL
RUNKUT(XA,Y,XB,NEQN,TOLA,TOLR,HSTART,WORK, &
IMETH,IERROR,ICOM,FUNC)
and must be called from a main "calling" program.
1.2. VARIABLES PASSED TO RUNKUT.
The variables passed are explained below and are
listed in the order in which they appear in the
call statement. The superscripted numbers refer to
notes appearing at the end of this section.
XA (Real-Input) XA is the starting
point of the interval over which integration of
the ODE's is to be performed.(1)
Y (Real array of size NEQN - Input/Output)
On entry, Y must contain the initial values of the
Y-variables--that is Y must contain the solutions at the
point XA. On exit from RUNKUT, Y will contain the solutions
at the end of the specified interval--that is at XB.(2)
XB (Real-Input) XB specifies the end
of the interval over which integration is to be
performed--the point at which solutions to the
differential equations are desired.(1,2)
NEQN (Integer-Input) NEQN specifies
the number of coupled differential equations in
the system to be solved.
TOLA (Real-Input) TOLA is the
absolute error tolerance required on the
solution.(3)
TOLR (Real-Input) TOLR is the
relative error tolerance required on the
solution.(3)
HSTART (Real-Input/Output) On the
very first call to RUNKUT in a calling sequence,
HSTART must contain a non-zero value as an initial estimate
for the step-size.(1) The actual value used is relatively
unimportant as the program will optimize step-sizes to keep
the tolerances within specified limits and keep the number of
function evaluations to a minimum. It is left to the user's
discretion to choose a value that will not require excessive
readjustment by the program. On exit from RUNKUT, this
variable will contain the last step size used by the program.
It is recommended that for subsequent calls to RUNKUT, the
step size returned by the program be used as the value for
HSTART to solve the equations over the immediately following
interval.
WORK (Real array of minimum size NEQN*17 - Input)
This is a work area made available to RUNKUT by the calling
program.
IMETH (Integer-Input) Indicates to
the solver, the order of the method to be used.
IMETH can take on the following values:
1. - Fourth order Fehlberg method
2. - fifth order Verner method 3. - seventh
order Verner method
While solving a given set of differential
equations it is possible to change the order of
the method over sucessive intervals simply by
setting IMETH to the desired value and
resetting ICOM(1) to 0 on each subsequent call to RUNKUT in
which a change of order is required.(4)
IERROR (Integer-Output) IERROR is an
error status indicator returned by RUNKUT after
each call. It can have one of the following values:
0 - no error. 1 - the sum of
TOLA and TOLR is less than 100 times machine
epsilon. This can be caused by setting both TOLA and TOLR
to 0, or by specifying values that are too small for the
program to handle sucessfully. In this case TOLR is set
to a default value. 2 - the
problem is too stiff for the method to handle, or
there is a discontinuity in the function. It is
impossible to proceed in either case.(2) 3 -
both conditions 1 and 2 above have occured.
IERROR is set to 0 on every entry to RUNKUT.
ICOM (Integer array of size 4) ICOM
is a communications vector. Each array element is used
to pass a specific item of information (as defined below) to
and from the subroutine.
ICOM(1) (Input) ICOM(1) must be set
to 0 on the first call to RUNKUT. ICOM(1) is
internally set to a non-zero value by the program
after the first call. For subsequent calls to RUNKUT with the
same set of equations, ICOM(1) must be left at this non-zero
value. It is possible to change the order of the method used
(see also IMETH) over sucessive intervals. In this case
ICOM(1) must be reset to 0 each time the order of the method
is changed.(4)
ICOM(2) (Input) A flag that indicates
to the program whether to perform checks on the
specified values of TOLA and TOLR.
0 - no checking of the error tolerances.
1 - program checks if the sum of TOLA and TOLR is at least
100 times machine epsilon. If not, TOLR is set to a
default value and the error indicator IERROR is set to 1.
If ICOM(2) is set to 0, the program assumes
that it will be receiving non-zero values of at
least one of TOLA or TOLR. Failure to check
this before entry to RUNKUT could result in
fatal program errors. Also, as will be shown in an example,
decreasing the error tolerances does not always result in
improved solution accuracy.
ICOM(3) (Input) A flag that indicates
to the program whether to use the default
values of mimimum and maximum step size set by the
subroutine.(5)
0 - default values of the maximum and minimum
step size are used. 1 -
allows the user to set values for the maximum and minimum
permissible step size. These values are set by including
the following FORTRAN common block statement in the main
"calling" program:
COMMON /CONS/HMIN,HMAX
ICOM(4) (Output) Status flag
that indicates the presence of round-off
error in the solution. 0 - No round-off
error. 1 - Severe round-off error possible
in answer.
ICOM(4) is reset to 0 on entry to RUNKUT.
FUNC (Function name) FUNC is
replaced by the name of the subroutine that is
part of the main user routine in which the differential
equations are specified (See section on user supplied
routines). The subroutine name must be declared external
by the following FORTRAN statement:
EXTERNAL [SUBROUTINE NAME].
1.2.1. NOTES
1. XA could be greater than XB. In other words
it is possible to use RUNKUT to solve the
differential equations in a negative
X-direction given initial values at XA. In
that case the step-size will also be negative.
HSTART could be specified as a negative value, but that
is not essential as the program internally adjusts the
algebraic sign of the step-size.
2. If a discontinuity or extreme stiffness is
encountered, RUNKUT will return in XB, the
value of X closest to the point of
discontinuity or stiffness, at which the
solutions are still within the specified error bounds.
The array Y will then contain solutions at this value of
X.
3. If the sum of TOLA and TOLR is less than
100 times machine epsilon (for example if
neither were defined in the main program),
TOLR is set to a default value of 10(6)
times machine epsilon. TOLA is not set to any
default value. The program uses a combination of TOLA and
TOLR to determine the accuracies of the solutions, and
specifying TOLR to a certain value almost always gives
better results than specifying TOLA to the same value.
4. Changing the order of the method on every
subsequent call to RUNKUT with a given
system of equations can lead to excessive
computaion times since a large number of
constants are re-evaluated each time the order is
changed.
5. The default values used by the subroutine
are 10(-8) for HMIN, and for HMAX as
follows:
0.5 - if the fourth order method is used.
1.0 - if the fifth order method is used.
2.0 - if the seventh order method is used.
However, if any of these values are larger
than the absolute value of the interval
width, HMAX is set to the interval width.
It is sometimes necessary to specify HMAX
to a small value to keep the step sizes to relatively
small values that keep the runge-kutta methods within
regions of stability. The subroutine detects the onset of
problem "stiffness" or the occurence of a discontinuity
in the function by checking if the step size is smaller
than HMIN. If the equations are being solved in the
negative X-direction, HMIN and HMAX are the smallest and
largest permissible step sizes in the negative
X-direction respectively. In other words, HMIN and HMAX
always define absolute constraints on the step size.
6. Machine epsilon on the IBM PC is very much
a function of the compiler used (even with
an 8087 math coprocessor installed). For
example Microsoft's Fortran-77 compiler
(with 8087 support) generates an epsilon of about
10(-16), whereas IBM/Ryan-McFarland's Professional
Fortran-77 generates an epsilon of about 10(-20)
1.3. THE USER-SUPPLIED ROUTINES
Two routines must be supplied:
1. A main "calling" program that calls RUNKUT.
It must also define the initial values of Y
at XA, set TOLA and TOLR, set ICOM(1) to 0
on the first call to RUNKUT, define the
number of equations in the sytstem to be solved,
dimension a WORK array for RUNKUT, and check for the
error status as passed back from RUNKUT through IERROR.
See the examples at the end of this section.
2. A subroutine containing the system of
differential equations. The subroutine
name must be declared EXTERNAL in the main
calling program and its name passed to RUNKUT
through the parameter FUNC-see above. The subroutine is
defined using the follwing FORTRAN statement:
SUBROUTINE [FUNC] (X,Y,YPRIME,NEQN)
where the parameters have the following
definitions:
X (Real) Contains the
current value of X as set in RUNKUT. Do
not alter this value within [func]
Y (Real - array of size NEQN)
Contains the current solutions at X. Do not alter in
[func].
YPRIME (Real - array of size NEQN)
YPRIME is set in the subroutine to the differentials
dy/dx as follows:
YPRIME(1) = F/1/(X,Y(1),...,Y(NEQN))
. .
YPRIME(NEQN)=F/NEQN/(X,Y(1),...,Y(NEQN))
where F/i/ is the i'th differential
function.